Preliminary work to compare the spatial cimis values at the cimis locations to the station’s actual reported values.
This assumes that you already have a local copy of ssj-weather cloned. If not get it from https://github.com/ssj-delta-cu/ssj-weather
new_stations <- read.csv("data/spatialcimis/cimis_newdeltastations_wy2016.csv", stringsAsFactors = FALSE)
ns <- new_stations %>% select(Stn.Id, Date, cimis_eto_new=ETo..mm.)
ns$Date <- as.Date(ns$Date, "%m/%d/%Y")
ns$Stn.Id <- as.numeric(ns$Stn.Id)
join_df <- dplyr::inner_join(cimis_data, sp_cimis_data, by=c("Date"="date", "Station"="id"))
join_df2 <- dplyr::left_join(join_df, ns, by=c("Date", "Station"="Stn.Id"))
cimis_both_eto <- join_df2 %>% select(Station, Date, cimis_eto_api=DayEtoValue, spatial_cimis_eto, cimis_eto_new)
cimis_both_eto <- cimis_both_eto %>% mutate(cimis_eto=ifelse(is.na(cimis_eto_api), cimis_eto_new, cimis_eto_api))
cimis_both_eto <- cimis_both_eto %>% select(Station, Date, cimis_eto, spatial_cimis_eto)
cimis_both_eto_long <- gather(cimis_both_eto, 'Type', 'eto', 3:4)
ts <- ggplot(cimis_both_eto_long, aes(Date, eto, colour=Type, group=Type))+ geom_line(size=1.15, alpha=0.7)+scale_x_date()+ylab("ETO (mm/day)")+
theme_bw()+
scale_color_manual(labels=c("CIMIS ETo", "Spatial CIMIS ETo"), values=c("cimis_eto"="red3", "spatial_cimis_eto"="royalblue1"))+
theme(panel.border = element_rect(colour = "black", fill="NA", size=1),
panel.grid.major = element_blank(),
plot.title = element_text(hjust = 0.5),
panel.grid.minor = element_blank(),
legend.position="bottom", # position of legend or none
legend.direction="horizontal", # orientation of legend
legend.title= element_blank(), # no title for legend
legend.key.size = unit(2, "cm"),
axis.text.y = element_text(size=12),
axis.text.x = element_text(size=12),
axis.title = element_text(size=14))+
guides(colour = guide_legend(nrow = 1),
fill = guide_legend(keywidth = 1, keyheight = 1),
linetype=guide_legend(keywidth = 1, keyheight = 1))+
facet_wrap(~Station, ncol=3)+theme(strip.background = element_blank(),
strip.placement = "outside",
strip.text.x = element_text(size = 16))
ts
station_timeseries <- function(station_num){
cimis_both_eto_long_single <- cimis_both_eto_long %>% filter(Station == station_num)
ts <- ggplot(cimis_both_eto_long_single, aes(Date, eto, colour=Type, group=Type))+ geom_line(size=1.1)+scale_x_date()+ylab("ETO (mm/day)")+
theme_bw()+
scale_color_manual(labels=c("CIMIS ETo", "Spatial CIMIS ETo"), values=c("cimis_eto"="red3", "spatial_cimis_eto"="royalblue1"))+
theme(panel.border = element_rect(colour = "black", fill="NA", size=1),
panel.grid.major = element_blank(),
plot.title = element_text(hjust = 0.5),
panel.grid.minor = element_blank(),
legend.position="none", # position of legend or none
legend.direction="horizontal", # orientation of legend
legend.title= element_blank(), # no title for legend
legend.key.size = unit(2, "cm"),
axis.text.y = element_text(size=12),
axis.text.x = element_text(size=12),
axis.title = element_text(size=14))+
guides(colour = guide_legend(nrow = 1),
fill = guide_legend(keywidth = 1, keyheight = 1),
linetype=guide_legend(keywidth = 1, keyheight = 1))
ts
}
stations <- unique(cimis_both_eto_long$Station)
for(s in 1:length(stations)){
print(stations[s])
p <- station_timeseries(stations[s])
print(p)
}
## [1] 47
## [1] 212
## [1] 167
## [1] 140
## [1] 242
## [1] 243
## [1] 247
## [1] 248
## [1] 249
sc <- ggplot(cimis_both_eto, aes(cimis_eto, spatial_cimis_eto, colour=factor(Station), group=factor(Station)))+geom_point(alpha=0.1)+geom_smooth(method='lm', se = FALSE)
sc<-ggplotly(sc)%>%
layout(autosize = F, width = 800, height = 600)
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
## Warning: Specifying width/height in layout() is now deprecated.
## Please specify in ggplotly() or plot_ly()
sc